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We consider the real time dynamics of an initially localized distinguishable impurity injected into 
the ground state of the Lieb-Liniger model. Focusing on the case where integrability is preserved, 
we numerically compute the time evolution of the impurity density operator in regimes far from 
analytically tractable limits. We find that the injected impurity undergoes a stuttering motion 
as it moves and expands. For an initially stationary impurity, the interaction-driven formation 
of a quasibound state with a hole in the background gas leads to arrested expansion - a period of 
quasistationary behavior. When the impurity is injected with a hnite center of mass momentum, the 
impurity moves through the background gas in a snaking manner, arising from a quantum Newton’s 
cradle-like scenario where momentum is exchanged back-and-forth between the impurity and the 
background gas. 


Introduction .— With recent experimental advances in 
the field of cold atomic gases, the physics of the one¬ 
dimensional Bose gas is receiving an increasing amount of 
attention.These systems, in which one has unprece¬ 
dented isolation from the environment and fine control of 
interparticle interactions, are excellent tools for examin¬ 
ing novel phenomena arising from strong correlations. 
One such phenomenon which has piqued both 

theoretical®^^. 16-21 experimentaR®’^^ curiosity is the 

expansion dynamics of a gas of cold atoms. A number 
of surprising and interesting effects have been observed, 
of which arrested expansion (or self trapping) is of 
particular relevance to this work. A gas (bosons^® or 
fermions^®’^^) is released from a confining potential and 
allowed to expand on the lattice. Under this time evo¬ 
lution ‘bimodal’ expansion is observed: the sparse outer 
regions of the cloud rapidly expand whilst the dense cen¬ 
tral region spreads only very slowly. This can be par¬ 
tially understood by considering the limit of strong in¬ 
teractions: doubly occupied sites are high energy config¬ 
urations which, thanks to the lattice imposing a finite 
bandwidth and energy conservation, cannot release their 
energy to the rest of the system and decay.®® 

With these recent experimental advances^® has also 
come the ability to examine systems in which there is 
a large imbalance between two species;®’®’’®’^®^^® a natu¬ 
ral starting point for the study of impurity physics. This 
gives insight in to a diverse range of problems,^® from the 
physics of polarons®'’’®® to the x-ray edge singularity®®’®® 
and the orthogonality catastrophe.®® The physics of im¬ 
purities also plays an important role in the calculation of 
edge exponents in dynamical correlation functions®^ and 
in understanding the nonequilibrium dynamics following 
a local quantum quench.®®’®®’®® 

The experimental study of the out-of-equilibrium dy¬ 
namics of a single impurity in the one-dimensional Bose 
gas has revealed rather rich physics: from how an impu¬ 


rity spreads when accelerated through a Tonks-Girardeau 
gas,® to how interactions effect oscillations in the size of a 
trapped out-of-equilibrium impurity.® Numerous theoret¬ 
ical investigations have addressed the Tonks-Girardeau 
regime: from a static point impurity®'’’®® to a completely 
delocalized (e.g., plane wave) impurity. ®®“^^ Away from 
the Tonks-Girardeau limit, theoretical study of the con¬ 
tinuum problem is challenging and results have focused 
on lattice models, such as the Bose-Hubbard model. ®’^® 


In this letter we consider the out-of-equilibrium dy¬ 
namics of an initially localized impurity in the Lieb- 
Liniger model. Using a combination of exact analytical 
results and numerical computations, we show that an im¬ 
purity injected into the ground state of the Lieb-Liniger 
model undergoes a stuttering sequence of rapid move¬ 
ment/expansion followed by arrested expansion. For 
an initially stationary impurity, this is caused by the 
interaction-driven out-of-equilibrium formation of a qua¬ 
sibound state of the impurity with a hole in the back¬ 
ground gas. This quasibound state is robust under time- 
evolution for long periods of time. For an impurity with a 
finite initial center of mass (GOM) momentum, the stut¬ 
tering sequence results in the impurity “snaking” through 
the background gas; the impurity exchanges momentum 
back-and-forth with the background gas through a quan¬ 
tum Newton’s cradle-like mechanism.® The results we 
present are relevant to experiments (see, e.g. Ref. [6]) 
and should be observed under reasonable conditions. 


The two-component Lieh-Liniger model .— We consider 
two species of delta function interacting bosons confined 
to a ring of length L. The Hamiltonian of the two- 


2 


component Lieb-Liniger model (TCLLM) is given by 



where herein we set h = 2m = 1, c is the interaction 
parameter, and the boson operators obey the canonical 
commutation relations (x), (j/)] = Sj^iS(x — y) with 

j,l = 1,2 denoting the species. As in the case of the one- 
component Lieb-Liniger model,the generalization to 
multiple particle species remains integrable provided all 
species interact identically. 

The TCLLM can be solved by the Bethe Ansatz,^®’^° 
giving access to some of its basic physical properties (see, 
e.g., [51-53] and references therein). An 7V-particle eigen¬ 
state containing 7Vi particles of species 1 is characterized 
by a set of N momenta {gjjv = {?!, • ■ ■ ,g 7 v} and a set 
of Ni species rapidities = {Ai,..., Aat^}. These 

momenta and rapidities satisfy the nested Bethe ansatz 
equations 
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where j = 1,... ,N and k = 1,... ,Ni. The eigenstate 
|{9}Af; {AIatJ has energy Eg = J2j q'j momentum 
Kq = X) j qj ■ 

The initial state. — We study the dynamics of an im¬ 
purity starting from the state 


l'i'(Q)) 




dx 


aQ 




( 4 ) 


where [fl) is the N 2 = N — 1 particle ground state of 
the one-component Lieb-Liniger model, Q is the COM 
momentum of the impurity and Af normalizes the state. 
The study of such a state is partially motivated by the 
experiments performed in Refs. 6 and 9, which study the 
dynamics of an impurity in a background gas. 

The initially localized impurity of Ref. 6 is prepared by 
illuminating a trapped one-component Bose gas with a 
radio-frequency pulse; this causes transitions between the 
|F, mi?) = |1, —1) hyperfine state of the trapped gas and 
the |1,0) state (the impurity). Due to the magnetic trap, 
transitions occur only within a spatially localized region, 
the thinness of which is Fourier-limited by the pulse du¬ 
ration. The resulting impurity contains up to three par¬ 
ticles and is accelerated through the gas by gravity, as 
the 11,0) state does not experience the magnetic trap. 

On the other hand, the impurity in Ref. 9 is prepared 
by first tuning the interspecies interaction to zero and 
then using a species-dependent trap and light blade to 


shape the impurity. Following this preparation, the inter¬ 
species interaction is turned on and the impurity released 
from the trap/light blade and its expansion studied. 

To distill the intrinsic dynamics of the impurity, our 
scenario varies slightly from experiments®’®: we study 
an impurity injected into a constant density background 
gas in the absence of an external potential (such as a 
magnetic trap and gravity). Similar approximations have 
been applied in the well-studied yrast states. 

Time evolution protocol. — Our aim is to compute 
the impurity density profile when the initial state (4) is 
time-evolved according to the Hamiltonian (1) pi{x,t) = 
{d>{Q)\e^^*d>\{x)'^i{x)e~'‘^*\d>{Q)). This is a nontrivial 
problem as the initial state (4) is not an eigenstate of the 
Hamiltonian. We use the integrability of the TCLLM to 
numerically evaluate the density profile using recently de¬ 
rived results for matrix elements of local operators.®® Due 
to a dearth of results for matrix elements in the TCLLM, 
we are restricted to studying the density of the impurity 
and we cannot examine the background gas.®® 

The essential idea is the following: we insert complete 
sets of eigenstates between each time evolution operator 
and the initial state in pi{x, t). By orthogonality, we sum 
over the Bethe states with Ni = \ and N 1 +N 2 = N. The 
momenta and rapidities characterizing these states sat¬ 
isfy the nested Bethe ansatz equations (2,3). The density 
profile of the impurity will then be given by 

Pi{x,t) = ^ (4'(Q)|{p}; A)({p};A|4'l(0)«'i(0)|{fc};^) 
{p}A 

X {{k}- (5) 

where {g} = The overlap of the initial state with 

a Bethe state can be expressed as A^dfe}; p,|4>((5)) = 

/(fdx ({fc};/r|4'l(0)|H), where 

Kq is the momentum of the ground state |H). So, in or¬ 
der to compute (5) we require two ingredients: the matrix 
elements of the creation operator 'I'|(0) and the density 
operator 4'|(0)4'j^(0) on the Bethe states. These ma¬ 
trix elements have been derived from the algebraic Bethe 
ansatz.®®’®® Required results are summarized in the Sup¬ 
plemental Materials.®® 

Readers interested in our scheme for numerically eval¬ 
uating the expansion (5) can consult Ref. [62]. An impor¬ 
tant point to note is that the expansion (5) contains an 
infinite number of terms. We truncate the Hilbert space 
by selecting the Bethe states which have the largest over¬ 
laps with the initial state (4). To quantify the truncation 
error, we compute the saturation of the sum rule 
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and we present numerical values for this with our results. 
We are limited to small numbers of particles ^ 10 and 
we have to keep ^ 10® — 10® states to saturate the sum 
rule to 2 decimal places. 
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FIG. 1. Time evolution of the impurity density of the initial 
state (4) with Q — 0, xg = I//2 and Og = 1.125 on the 
L = 40 ring for a system of = 8 particles with interaction 
parameter c = 10. The Hilbert space is truncated to 25150 
states, leading to the sum rule (6) = 0.9858. (Inset) Time- 
evolution of the maximum of the density pi(xo,t). Constant 
time cuts can be found in the Supplemental Materials.®^ 



FIG. 2. The dynamics of the impurity density (5) from the 
initial state (4) with Og = 1.125 and Q = 407r/L. We use 
21,507 states to study a system of = 8 bosons on the length 
L = 40 ring with interaction parameter c = 10, resulting in 
the sum rule (6) = 0.981. Plots of constant time cuts are 
presented in the Supplemental Materials. 


The noninteracting limit .— In the noninteracting 
limit, the time evolution of the initial state (4) is a 
single particle problem. The time-dependent density 
profile can be calculated exactly (we take L —>■ oo): 

Pi{x,t)c=Q = ^ exp(-)/\/ao + The non¬ 
interacting density profile remains Gaussian at all times, 
with a time-dependent width and amplitude. 

Arrested expansion: Q = 0 .— In Fig. 1 we present re¬ 
sults for the time evolution of the impurity density pro¬ 
file (5) for N = 8 bosons on the circumference L = 40 
ring starting from the initial state (4) with Xq = T/2, 
Uq = 1.125 and interaction parameter c = 10. We mea¬ 
sure time in units of tp = l/Ep where Ep = is 

the Fermi energy in the c ^ oo limit. The Hilbert space 
is truncated to 25150 states, leading to the sum rule sat¬ 
uration 0.9858 (i.e., to 1.4%). Upon time evolution the 
wave packet spreads, maintaining its Gaussian shape as 
in the noninteracting case. However, at time t ^ 2tp the 
wave packet stops spreading and only undergoes small 
amplitude breathing oscillations. This arrested expan¬ 
sion is an example of prethermalization.^^^'^^ The system 
relaxes in a two-step process, first approaching a quasi¬ 
stationary non-equilibrium state (the arrested expansion) 
before subsequent equilibration. Two-step relaxation has 
been observed in the one-dimensional Bose gas following 
a global quantum quench. 

We can qualitatively reproduce aspects of this behavior 
with a mean field (MF) decoupling of the interaction term 

4'](a;)4'|(x)4'j (x)4'^ (x) « pj(x, t)4'J(x)'I'; (x) +j O 1. 

(7) 

At a MF level the impurity profile is a time-dependent re¬ 
pulsive one-body potential for the background gas. The 


region under the impurity then excludes particles in the 
background gas, resulting in the formation of a ‘hole’. 
This hole in the background gas acts as a confining (at¬ 
tractive) one body potential for the impurity in MF and 
the two form a quasibound particle-hole pair,^^ much like 
an exciton in the electron gas (see, e.g.. Ref. [75]). This 
is different to the self-trapping scenario on the lattice: 
there the ‘doublons’ are stable as the large interaction 
energy cannot be converted into kinetic energy due to 
particle number conservation and the finite bandwidth. 
In the continuum, dynamical arrest is driven by the for¬ 
mation of the impurity-hole quasibound state and is not 
observed for an indistinguishable impurity.®^ 

At later times (t > 7tp), the impurity eventually 
broadens. This broadening occurs in a sequence of expan¬ 
sion/arrested expansion steps, whilst the impurity under¬ 
goes small amplitude breathing oscillations.®^ The slow 
decay of the density at later times may be related to the 
subdiffusive equilibrium behavior reported in Ref. [35]. 
However, finite-size effects and our choice of observable 
obscure the characteristic logarithmic decay of subdiffu- 
sion. 

The snaking impurity: Q ^ Q .— Finally, we consider 
the time evolution of the initial state (4) with nonzero 
GOM momentum Q. Our prescription for computing 
the time evolution is identical to the Q = 0 case; in 
Fig. 2 we present results for the impurity density profile 
for the same set of parameters as in Fig. 1 with Q = n. 
We see rather surprising behavior: the impurity moves 
in a snaking manner, repeatedly moving and expanding 
before becoming approximately stationary with arrested 
expansion. To quantify the nonuniform motion of the 
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impurity further, we define the COM coordinate X{t) as 


X{t) = — arctan 

ZTT 


d0 sm6pi{6,t) 
Jq^ d0 cos 0pi(0,t) 


( 8 ) 


where 0 = 2'kxIL. We plot the COM coordinate in Fig. 3 
for a number of interaction strengths; X{t) shows regions 
of rapid movement, followed by (approximately) station¬ 
ary plateaux. Only at t < tp/3 does the COM move as 
in the noninteracting case: X(t)c=o = -^(0) — 2Qt. The 
sharpness of the plateaux and transient regions are gov¬ 
erned by the interplay between the delocalization of the 
impurity and its interactions with the background gas.^® 



Time t/tp 

FIG. 3. The time evolution of the center of mass X{t) (8) 
for: (points) the initial state (4) with Aq = 1.125, Q — tv for 
N = 8 particles on the length L = 40 ring with interaction 
parameter c = 5,10, 20; (line) a noninteracting point particle 
with mass m = 1/2 and velocity Qlm. Inset: velocity of the 
center of mass V{t) — AX{t)/At [cf. Ref. 39). 

We have the following picture for the behavior shown 
in Fig. 2: (i) The impurity moves to the left, scattering 
particles in the background gas and creating excitations 
with finite momentum, (ii) The impurity continues to 
scatter with the background until it imparts most (or 
all) of its COM momentum, (iii) The excitations in the 
background gas propagate around the ring and then col¬ 
lide once more with the impurity, (iv) The impurity gains 
COM momentum and the process repeats. In support of 
this picture is the behavior of the COM position plateau 
with system size L: the time for leaving the plateau Tp 
is (approximately) linearly dependent on L. Tp is also 
related to the initial momentum Q of the impurity; for 
large Q, Tp ~ 1/Q (an excitation with momentum Q has 
velocity ~ Qjm). This Q-dependence reflects the mo¬ 
mentum imparted by the impurity to excitations in the 
background gas, which then propagate around the ring.^^ 
This process can be thought of in terms of a quantum 
Newton’s cradle^ on a ring, with the impurity exchang¬ 
ing momentum back-and-forth with the background gas, 
resulting in the snaking motion shown in Fig. 2. In the 
Supplemental Materials we show that this behavior is 
not realized on the lattice when we perform the MF de¬ 


coupling (7) for the same set of parameters that capture 
some aspects of the Q = 0 behavior. 

It is interesting to consider removing periodic bound¬ 
ary conditions: excitations produced by injecting the im¬ 
purity will propagate towards the boundary and subse¬ 
quently reflect, returning to once again scatter the im¬ 
purity. This reflection of the excitations means that we 
expect the COM to snake back-and-forth about xq rather 
than around the ring. In the presence of a harmonic trap, 
the COM will travel in a snaking motion due to both the 
trap and collisions with the background excitations. 

A question that has recently attracted attention 
is whether an injected impurity has finite momen¬ 
tum in the t —>■ oo limit (see, e.g.. Refs. [39- 
44]). To address this, we compute the momen¬ 
tum of the impurity in the diagonal ensemble 

(de)78 Koe = E{fc};M('i^(Q)IW;M)({fc};Ml^(Q)) x 
J 2 p{{k}-, 0 \p'd>lpdfij{k}-p), where T-pp 
1 /L J da:e“*P’^'I'i(a:). Doing so, we find Kjjp « —0.022 
(for N = 4 particles on the length L = 40 ring), in 
keeping with general expectations from the study of the 
delocalized impurity in the Tonks-Girardeau limit. 

We have also examined the density of the impurity in 
the DE to ascertain whether translational symmetry is 
restored in the long-time limit. Generically, we find that 
translational symmetry is not restored in the finite-size 
system due to a symmetry of the Bethe states under a 
change in sign of all the momenta and rapidities. 

Conclusion .— In this letter, we consider the nonequi¬ 
librium time evolution of a single localized impurity (4) 
injected into the ground state of the Lieb-Liniger model. 
In both the case of zero and finite COM momentum, we 
observe a ‘stuttering’ behavior in the motion. In the first 
case (see Fig. I), this quantum stutter manifests in the ar¬ 
rested expansion of the impurity (in the absence of a lat¬ 
tice). This arises from the out-of-equilibrium formation 
of quasibound impurity-hole pairs which are stable for 
extended periods of time. This interaction-driven effect 
can be qualitatively captured by the MF decoupling (7): 
the impurity repels the background gas, leading to the 
formation of a hole which acts as a confining potential for 
the impurity. Eventually the impurity broadens in a se¬ 
quence of rapid expansions and quasistationary periods, 
all the while undergoing small amplitude breathing oscil¬ 
lations. This stuttering motion and the quasibound state 
formation highlights the importance of distinguishability, 
as this mechanism does not exist for an impurity of the 
same species as the background gas.^® 

In contrast, when the impurity is injected with a finite 
COM momentum, the quantum stutter is clearly seen in 
the motion of the impurity, which snakes through the 
background gas, see Figs. 2 & 3. We can picture this 
as a quantum Newton’s cradle^ on the ring: the injected 
impurity scatters particles in the background gas until it 
loses most of its COM momentum. These scattered exci¬ 
tations then propagate around the ring and subsequently 
collide with the impurity, causing it to move once again. 
This process repeats, leading to the stuttering, snaking 
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motion of the COM. Quantum flutter, the exchange of 
momentum back-and-forth between a delocalized impu¬ 
rity and the background gas, has been studied in the 
Tonks-Girardeau regime.The momentum of the im¬ 
purity in the long-time limit was computed by means of 
the DE and found to be small, but non-zero. 

Our results are of direct relevance to experiments 
in cold atomic gases and the observed physics should 
not be reliant upon the integrability of the model 
(see, e.g.. Refs. [40] and [61]) and should survive finite 
temperature.The discussed results may also be use¬ 
ful in elucidating the properties of the TCLLM at finite 
temperature, where it is likely that impurity-like solitons 
arise.Finally, this work provides a nontrivial check 


and validation of cutting-edge theoretical results for the 
matrix elements of the TCLLM in the extreme imbalance 
limit.53.60 
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Supplemental Materials for “Motion of a distinguishable impurity in the Bose gas: 
Arrested expansion without a lattice and impurity snaking” 


A. Summary of required matrix elements for the 
two-component Lieb-Liniger model 

Here we summarize known results®^ for matrix ele¬ 
ments of local operators in the two-component Lieb- 


Liniger model. For technical reasons, the known matrix 
elements of Ref. [53] are restricted to the case with a 
single impurity boson (Ni = 1) and to the local oper¬ 
ators 'I'i(O), T'|(0)T'i(0) and 'I'|(0)'I'2(0). As a conse- 
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quence of this, we are not able to ‘image’ the background 
gas 4'2(0)'I'2(0)- Very recently®^ there have been new 
(and more general) results for matrix elements in the 
two-component Lieb-Liniger model, but we have yet im¬ 
plement them. 


1. Normalization conventions 

We focus on Wparticle eigenstates containing a single 
Ni boson. For clarity and ease of comparison, here we 
work with the conventions of Ref. [83] and we define the 
(non-normalized) eigenstates |{q}Ar;/x)). The eigenstates 
in the body of the text are recovered by normalization 
\{q}N;iJ-) = |{g}jv;M))/ll{9}Af;Mll> where norms of the 
eigenstates are given by 


II{9}jv;mIP = fJ.)) = c detf72, 

J2 = 


Jqq Jqfi 
Jfj.q Jv-n 


Here is the Jacobian of the nested Bethe ansatz equa¬ 
tions [see Eqs. (2,3)[ given by the matrix elements 


{Jqq)jl — Sjl 


N 




- qi), 

“ ‘P2{kj — ^), 

N 


m—1 


where we define the scattering phase ipn{u) = 
2cn/(v?v?' + c^) and j,l = 1,..., N. 

We take the one-component Lieb-Liniger eigenstates 
[{pIat}) to have their conventional normalization^®’®^^®® 


IIWjvIP = n [(Pi -Pi) 

3<l 

N 


^ +C^ 


det Ji, 




PliPj - Pm) 


- PliPj -Pl)- 


2. Matrix elements (({p}]v|J'i(0)|{fc}jv; A)) 

For two states with no coinciding momenta, the matrix 
element of the impurity annihilation operator takes the 
determinant form 

^ n^>j{h-k,+ic) -^c 

Ul>TniPl - Pm + ic) nj(A - kj - ic/2) 


Here the (TV — 1) x [N — 1) matrix IM has elements Y\jk = 
Mjk — with 


M,k = t{pk - k,)h2{\ - 


X\Z.=ihi{km - kj) 

I ,/7 \U (U \\T\.m=l^t{kj — Pm) 

Pk)h 2 \kj A) ^ 

Y[m=l hl{kj — km) 

where we’ve defined the functions hn{u) = u + icjn and 
t(u) = —c/[u(u -\- ic)]. 


3. Matrix elements (({p}jv; (0)J'i(0)|{fc}jv; A)) 


For two states with no coinciding momenta, the matrix 
elements of the impurity density operator are 


{{{p}N;pmmim{k}N;x)) 
—i 
c 


= — (-l)^(^+i)/2TT.^-1- _TT. --_detV 

r J- i’• — hi — ir r) - — ' 


kj - ki-ic -I -pj -pi+ic 

J>1 J 


N 

n i’^i 


-Pm+ic 


n, (A - kj - ic/2){p-pj + ic/2) 


I,m—1 

where the (A^ + 1) x {N + 1) matrix V has elements 


ic 


ic 


^ji = - A -k - j M - -^y{k] - Pl) 

+ {pi-X-'^'){pi-p+t)i{Pi-k3) 


N 

n 


- fern + ic){pi -pm- ic) 
{Pl - km- ic){pi -Pm-\- ic) 


yAT+ij = 


N 

n 


Pm - Pj + ic 


km - Pj + ic ’ 


V 


j,N+l — Ij 


>N+1,N+1 


= 0 , 


with j,l = 1,...,TV. Additionally, the diago¬ 

nal elements of the impurity density operator fol¬ 
low from translational invariance of the eigenstates: 
({fc}^;A|vI/I(0)vI/i(0)|{fcU;A) = l/L. 


B. Dynamics of an indistinguishable impurity in 
the one-component Lieb-Liniger model 

Here we consider the initial state 

L 

|^2(Q)) = ^ ) 'I'J(x)|L!), (SI) 

0 

where jH) is the ground state of N 2 bosons of species 
2, e.g. the analogue of Eq. (4) with an indistinguish¬ 
able impurity. We wish to consider the time evolu¬ 
tion of the expectation value of the density operator 
















Gaussian shape with oscillations superimposed on top. 

We present the time evolution of the density opera¬ 
tor on the initial state (SI) for iV = 8 particles, <5 = 0, 
Oq = 1.125 and L = 2xq = 40 in Fig. SI. We see that the 
behavior of the density as it evolves in time is qualita¬ 
tively consistent with the noninteracting result (S2) . We 
see no evidence of a stalling of the spreading of the wave 
packet, which is not surprising as the mechanism which 
exists in the two-component case (hole formation in the 
background gas and subsequent trapping of the impurity 
in the hole) is not present when there is a single species 
of boson. 

C. Lattice mean field description 


FIG. SI. Time evolution of the density profile P 2 {x, t) for the 
initial state (SI) in the one-component Lieb-Liniger model for 
N = 8 particles on the circumference L = 40 ring with inter¬ 
action parameter c = 10. We used 12619 states to saturate 
the sum rule to 0.9937. 


2 {x) with this initial state. Our prescription 
for computing the time evolution is analogous to the 
two-component case: we expand the expectation value 
P 2 {x,t) = (4'2(Q,t)|4'J(a:)4'2(a:)|4'2(Q,t)) in terms of 
known matrix elements and overlaps between the Bethe 
states,^®’®^^®® to obtain an expansion similar to Eq. (5). 

For the case with Q = 0, we can gain some insight 
from examining the noninteracting limit c = 0. Working 
on the infinite system (L —> oo) with xq = 0, we find the 
density 


P2{x,t) = p + 


2'Ka^p 


e 0 ^ 


\/al+ 4t2 


-I- p\Ph^ 


K+4f®)i 


r cos 6»t - 


tx^ 


■ 4^2 


(S2) 


where p is the average density and with p and 29t = 
arctan(2t/ao). Thus we expect the wave packet to be of 


In an attempt to explain the dynamics of the initial 
state, Eq. (4) in the main text, we consider the following 
lattice Hamiltonian (we consider a lattice Hamiltonian 
for numerical convenience): 

.Hiatt = ~Jd ^ ^ + H.c^ -|- U ^ ) {n\ ~\- nf^ 

i i 

-J6E(^^+i+H-c.), (S3) 

I 

where = ajai is the number operator. This model 
is motivated by the two-component Lieb-Liniger model: 
we consider two species of bosons which have an on-site 
interaction only and the kinetic terms coincide in the 
continuum limit 

lim d{x)^d(x + oq) —>■ d'' {x)d{x) + {x)dxd{x) 

ao-tO 

+ {x)dld[x), 

'^Ad+i + d\^i-i const. + ao dx d)[x)dxd{x). 

i 

We choose the d bosons to play the role of the background 
gas (species 2 in the main body of the manuscript) and 
we consider U ^ Jd to reflect the strong coupling regime 
of the main text. 

The Heisenberg equations of motion for the boson bi- 
linears take the form 




- blb^^^ - -b iU 6l6j(n^ - n^^) -b (n^ - n’^^)blbj + 2{nf - nf)blb^ 


and similarly for djdj with d b. We take the expectation value of this expression and perform a time-dependent 
mean-field decoupling which preserves the U{1) symmetry for each of the species (cf. Eq. (7)) 

{{nf-nj)b\b^) (^{nf{t)) - (n^(<))) (&I(t)fcj(t)), {n'l{t)b\{t)b^{t)) {n\{t)){b\{t)b^{t)) + {d^bJ{t)){b^b\{t)), 

to arrive at the approximate equations of motion for the boson bilinears: 


= -iJb[{bl_ib^it)) - (blbj^dt)) + (bl+ib^it)) - {blbj_dt)) 
+2tU\{nfit)) - (n^(t)} + 2{n\{t)) - 2{n){t))\ {b\bAt)), 


(S4) 
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with similar for the d bosons. Our initial conditions are 
fixed by the initial state |'I'((5)); for the purposes of con¬ 
venience, we consider the ground state jO) to be the c = 0 
ground state (this is an approximation, alternatively one 
can view this situation as a combination of injecting the 
impurity and performing a quantum quench of the inter¬ 
action parameter), in which only the zero-mode is popu¬ 
lated. The initial conditions for the bilinears are: 


{dUj)o = P, 





(S5) 


where p = N/L is the density of the background gas and 
lAfp = J2x exp[—(x — xo)^/ao] is a normalization factor. 



FIG. S2. Time depedence of the (a) ni,-, (b) rid boson number 
operator expectation values from the mean-field equations of 
motion (S4) with U = 14.5 and Jd = Jb = 1. Initial condi¬ 
tions (S5) with ao = 2, Q = 0 and p = 0.2 for L = 40 sites 
were used. 


(a) 


(b) 


FIG. S3. Time depedence of the (a) rib; (b) rid boson number 
operator expectation values from the mean-field equations of 
motion (S4) with U = 14.5 and Jd = Jb = 1- Initial condi¬ 
tions (S5) with ao = 2, Q = 13 X 27r/L and p = 0.2 for L = 40 
sites were used. . 



We present results for the expectation values of the 
density operators in Figs. S2 and S3 for U = 14.5 and 
Jd = Jb = ^ with p = 0.2 and ag = 2 on the circumfer¬ 
ence L = 40 ring. Parameters were chosen in an attempt 
to qualitatively reproduce the Q = 0 continuum behav¬ 
ior: arrested expansion followed by eventual spreading 
of the impurity (cf. Fig. S4). In Fig. S2 we see approx¬ 
imately the required behavior for Q = 0: the impurity 
initially spreads, but for times t tp — 9tp expansion is 
arrested (the nature of the amplitude and fluctuations is 
clearly very different in the mean field lattice case com¬ 
pared to the continuum) before subsequently spreading. 
In the background gas. Fig. S2(b), we see that a region 
of depleted density (a ‘hole’) appears below the impurity, 
which remains despite multiple collisions with propagat- 
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ing wave packets (the red peaks crisscrossing the figure). 

In Fig. S3 we present similar data for the case with 
Q = 13 X (27r/L) (we move away from Q = tt as this is 
a special point in the lattice case). Surprisingly we see 
that the addition of finite center of mass momentum for 
the impurity has lead to a strengthening of the dynamical 
arrest in the mean field approximation, with a deep and 
more robust hole forming in the background gas. Clearly 
we see no evidence of the ‘snaking’ behavior observed 
in the continuum, despite multiple collisions with excita¬ 
tions in the background gas. This strongly suggests that 
the behavior observed in the continuum for Q 7 ^ 0 is be¬ 
yond mean field theory (and may differ dramatically to 
that observed on the lattice). 

D. Addition plots: Q = 0 long time and Q = 0, tt 
constant-time cnts 

Here we present additional data for the time-evolution 
of the initial state with Q — 0 and Q = t:. In Fig. S4(a) 
we present the time-evolution for up to time t ~ 20tF- 
As we saw in the main text, for times ^tp < t < Itp 
the impurity undergoes arrested expansion: it is approx¬ 
imately stationary, with only small amplitude breathing 
oscillations. Following the arrested expansion, there is 
a period of rapid expansion, followed by a shorter qua¬ 
sistationary period and then subsequent expansion. In 
Fig. S4(b) we present constant time cuts for short times 
t < Q (t < ii.btp) , which show the initial period of rapid 
expansion and subsequent arrested expansion. Figure S5 
presents similar time cuts for the initial state with Q = tt 
for c = 5,10, 20. 


E. Quantum Newton’s cradle on the ring 

1. Motion of the center of mass 

Here we provide support for our assertion that the mo¬ 
tion of the center of mass coordinate can be explained 
in terms of a “quantum Newton’s cradle” on the ring. 
The motion of the center of mass coordinate is presented 
in Fig. 3 of the manuscript, where we observe periods 
of rapid motion separated by approximately stationary 
plateaux (see also Figs. S5). We have the following ex¬ 
planation for the observed behavior: as the impurity 
moves through the background gas, it collides and excites 
the background gas, with the excitations predominantly 
moving in the same direction as the impurity. The impu¬ 
rity continues to collide with the background gas until it 
has imparted all (or most) of its center of mass momen¬ 
tum, and subsequently the center of mass coordinate is 
(approximately) stationary. The excitations in the back¬ 
ground gas propagate around the ring, until they once 
again reach the impurity and collide with it, imparting 
momentum and causing the center of mass of the impu¬ 
rity to once more move. This process then repeats. In 


support of this picture, in Fig. S 6 (a) we present the cen¬ 
ter of mass motion in three systems with different sizes 
and fixed particle density (this means the speed of sound 
in the system should be similar, up to small finite size 
effects), and we see that the time for which the center 
of mass is stationary is linearly dependent on the system 
size L. In Fig. S 6 (b) we present the center of mass mo¬ 
tion for four different initial momenta Q, and we see the 
length of the plateau is inversely proportional to Q at 
large Q (the velocity of an excitation with momentum Q 
is V = Qjm). Both of these results are consistent with 
the presented picture, where the deviation from the sta¬ 
tionary plateau is driven by finite momentum excitations 
propagating around the ring. 


2. Behavior with variation of the interaction strength c 

How the behavior of the ladder motion of the cen¬ 
ter of mass changes with the interaction strength reveals 
the competition between interactions in the system and 
spreading of the impurity. An intuitive picture to have 
in mind is that of a liquid, which becomes ‘stiffer’ with 
increasing interaction strength. At first blush, such a 
picture may seem to be inconsistent with the presented 
center of mass motion (see Fig. 3 of the main text), which 
appears to sharpen with weakened interactions. The 
observed behavior can be explained as follows. When 
the interaction strength c is weak, the impurity spreads 
quickly and has almost completely delocalized by the sec¬ 
ond plateau. As a consequence, the second plateau is rel¬ 
atively fiat and the transient region between the plateaux 
is broad. With strengthening interactions, the spreading 
of the impurity is hindered and the spreading has yet to 
finish by the time the second plateau is reached. The sec¬ 
ond plateau appears less stable for strong interactions as 
the impurity continues to slowly spread whilst approx¬ 
imately stationary, shifting the center of mass slightly. 
We give evidence for this picture in the remainder of this 
section. 

The behavior of the center of mass motion can be seen 
in Fig. 3 of the main text for c = 5,10, 20 (see also Fig. S 6 
for c = 10). We will first address the behavior of the 
transient regions and then the stability of the plateaux 
upon varying the interaction strength, showing that it is 
consistent with the intuitive picture of interactions creat¬ 
ing a stiffer background gas. In the first transient region 
t < tp 12 the momentum of the impurity is imparted to 
the fluid: this happens more quickly for stiffer (c larger) 
fluids. In the second transient region, the impurity is 
accelerated by collisions with the excitations of the back¬ 
ground gas. The impurity immersed in the stiffer fluid is 
accelerated over a shorter period of time, and imparts its 
momentum back to the gas quicker, consistent with the 
intuitive picture. As a consequence, the transient region 
reduces in temporal extent with increasing interaction 
strength. 

Next, we consider the behavior of the second plateaux. 
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FIG. S4. Time evolution of the impurity density of the initial state |4'(Q)) [defined in Eq.(4)] with Q — 0, xo = Ljl and 
Oq = 1.125 on the L = 40 ring for a system of = 8 particles with interaction parameter c = 10. The Hilbert space is truncated 
to 25150 states, leading to the sum rule [Eq. (6)] = 0.9858. (a) The full time-evolution, showing the initial rapid expansion, 
followed by a period of arrested expansion and subsequent spreading/quasi-stationary periods. Inset is the time-evolution of the 
density at the midpoint, (b) Constant-time cuts at short times, showing in detail the rapid initial expansion and the period of 
arrested expansion. At intermediate times 5 < t < 18 (2tF TIf) the impurity is approximately stationary and Gaussian 

in shape. 


It is useful to consider Fig. S5; we see that the spreading 
of the impurity is suppressed with increasing interaction 
strength (this is particularly apparent in the second and 
fourth rows). We also see that the impurity has almost 
completely delocalized around the ring when c = 5. This 
is important, as the local fluctuations in the density are 
proportional to the local derivative of the density, and 
hence fluctuations are suppressed with increasing delo¬ 
calization of the impurity, which improves the stability 
of the plateaux. Increasing the momentum of the impu¬ 
rity also increases the rate at which the impurity delo¬ 
calizes, and the plateaux are flatter, see Fig. S7 [see also 
Fig. S6(b)). Furthermore, Fig. S5 also shows that the 
slow drifting of the second c = 10, 20 plateau corresponds 
to small changes in the shape of the impurity, correspond¬ 
ing to a transfer of weight leftwards, see Fig. S8. This 
transfer of weight, due to the slight asymmetric spreading 
of the impurity, causes the drifting of the second plateau. 

Delocalization of the impurity does not only increase 
the stability of the plateaux, it also leads to a smear¬ 
ing of the transient region, see Fig. 3 of the main text 
and Fig. S7(a). This is also consistent with our intuitive 
picture: excitations in the background gas now scatter 
on a increasingly extended object. Some excitations pass 
through the impurity, some scatter on the right or left 
of the impurity; the transient motion becomes more un¬ 
clear and the transition between plateaux broadens [e.g., 
the peak velocity of the COM is reduced, see the inset of 


Fig. 3 of the main text and Fig. S7(a)). 

F. The diagonal ensemble 

1. The impurity density in the diagonal ensemble 

To ascertain whether the impurity density becomes 
translationally invariant at long times after the impurity 
is injected, we compute the density profile in the diagonal 
ensemble 

Pi{x)de = Y. E 

{k}-n {p};A 

x(M;A|vI/t(0)vI/,(0)|{fc};Ar)({fc};Ar|4'(g)), 

(S6) 

which follows from a stationary phase argument in the 
long time limit. Herein we assume the diagonal ensem¬ 
ble coincides with the long time limit. Representative 
results for = 8 particles on the length L = 40 ring 
with interaction parameter c = 10 are shown in Fig. S9. 
There, we see that the diagonal ensemble result for the 
density with initial Q = 0 is not translationally invari¬ 
ant, whilst for Q = n the density profile appears much 
closer to constant. This provides strong evidence that 
for sufficiently large Q the impurity is almost completely 
delocalized around the ring in the long-time limit. 
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FIG. S5. Real time slices of the time-evolution of the impurity density for the three values of the interaction parameter 
c = 5,10, 20 from the initial state with ag = 1.125 and Q — n. 


One important question raised by Fig. S9 is whether 
there is a sharp transition or a smooth reduction in the 
extent to which translational symmetry is broken with 
increasing Q. In Fig. SlO(a) we present the impurity 
density in the diagonal ensemble for a number of ini¬ 
tial momenta Q and = 4 particles, and we show that 
the severity of the translational symmetry breaking is 
smoothly reduced as a Gaussian in the momentum of the 
initial state in Fig. SlO(b). Formally, this means that 
translational invariance is only recovered in the Q —1 oo 
limit. However, in a practical sense, translational symme¬ 
try is restored for sufficiently large Q for a finite precision 
measurement. 

The origin of this Gaussian scaling can easily 
be explained by a degeneracy in our system: for 
each Bethe state |{fci, ^ 2 ,..., ^at}; ir) there exists a 
state with the same energy and opposite momentum 
|{—fci, —fc 2 ,..., —fciv}; —fi). The diagonal ensemble (S6) 
now contains two types of terms: diagonal matrix ele¬ 
ments which sum to 1 /L and off-diagonal terms 5pi (x) de 


which break translational invariance: 

5pi{x)de = 5] (4'(Q)|{fc};M)({-fc};-M|'k(Q))e2*^''" 

+ (S7) 

Here the ellipses denote other terms arising from other 
(possible) degeneracies. The overlap between the initial 
state and the Bethe states (4>((5)|{/c};/r) is weighted by 
a Gaussian factor cx exp(—Q^). This Gaussian term will 
be present in any terms which break translational invari¬ 
ance, and hence the extent to which translational sym¬ 
metry is broken is smoothly suppressed with increasing 
Q, as observed in Fig. SlO(b). 

When considering the center of mass coordinate mo¬ 
tion, we do not see a qualitative difference for cases in 
which the diagonal ensemble result is (almost) transla- 
tionally invariant and those in which the translational 
invariance is more strongly broken, see Fig. S6(b). 




































13 


(a) (b) 




Time t/tp 


Time t/tp 


FIG. S6. Motion of the centre of mass coordinate X{t) [as defined in Eq. (8)] for the initial state [^(Q)) [defined in Eq. (4)[ 
with xo = L12 and Q = tt for (a) a number of system sizes L with fixed total particle density N/L and interaction strength 
c = 10; (b) a number of different initial center of mass momenta Q for L = 30, N = 6 and c = 10. 
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FIG. S7. (a) Motion of the center of mass of an impurity with Uq = 1.125 and Q = tt, 37r/2 for N = S particles on the length 
L = 40 ring with c = 10. (b) The impurity density on the (b) first and (c) second plateau for Q = 37r/2 (c/. the second column 
of Fig. S5 for Q = tt). 


2. The momentum of the impurity in the diagonal ensemble diagonal ensemble result is 


Having computed the diagonal ensemble result for the 
density of the impurity, we now turn our attention to 
computation of the momentum of the impurity. We de¬ 
fine the momentum of the impurity as 


K{t) =Y,P (S8) 

P 


where = r / is the momentum space 

annihilation operator for a boson of species 1. In the 
t oo limit we assume that this is given by the diago¬ 
nal ensemble and as ifI p^i.p conserves momentum, the 


Kde = EE p{^m{k};T){{khT\^m 

{fe};p P 

(S9) 

Fourier transforming to real space operators, and insert¬ 
ing the resolution of identity over one-component Lieb- 
Liniger eigenstates, we find 


Kde 


= H E - *■.) 

ffc};!* {«} 


{nQ)\{k};T){{k};p\^H0)\U}) ,(S10) 


which is expressed in terms of known matrix elements. 
Using the previously discussed symmetry of the Bethe 
states (and a similar property for the one-component 
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FIG. S8. The difference /S.pi{x,t\,t 2 ) = p\{x,t\) — pi{x,t 2 ) 
in the density at times t/tp = 4.34 and t/tp = 3.75 for c = 
20, see Fig. S5. We normalize to the maximum value of the 
density at time t/tp = 4.34. A clear transfer of weight, from 
the right hand side of the impurity wave packet to the left 
with increasing time is seen, resulting in the drifting of the 
center of mass plateau shown in Fig. 3 of the main text. 


states) and the properties of the matrix elements of the 
creation operators, we can write the momentum of the 
impurity in the diagonal ensemble as 


Kde “o(Q+^fc)^ _ e “o(Q 

{k}\fl {q} 


X 


(fl|4^(0)|{fc};M)({fc};A.|^nO)IM) 


(Sll) 


for the one-component ground state |n) containing an 
odd number of particles. We immediately see that for 
Q = 0 the momentum of the impurity remains at zero at 
all times (consistent with the Q = 0 time-evolution). 


FIG. S9. Results for the impurity density in the diagonal 
ensemble (S6) for the initial state |4'(Q)) [defined in Eq. (4)] 
with momentum Q = 0, tt and Uq = 1.125 for N = 8 particles 
on the length L = 40 ring with interaction parameter c = 10. 
The largest sum rule saturation requires 69532 states in each 
sum of Eq. (S6). Positions of the peaks are at xq and xo+L/2. 


(a) 



(b) 



FIG. SIO. (a) Deviation of expectation value of the impu¬ 
rity density operator in the diagonal ensemble (S6) from the 
translational invariant case for the initial state with momen¬ 
tum Q and xo = L/A. Results are for N = A particles on the 
ring of length L = 40 with interaction parameter c = 10; sum 
rules are saturated to at least 0.997 and we see now change in 
the extent of translational symmetry breaking with increas¬ 
ing saturation of the sum rule, (b) Maximum value of the 
deviation at a; = a:o as a function of Q with a Gaussian fit. 
Note that the peaks in (a) are shifted compared to Fig. S9 
solely due the change in xq. 

















